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INTRODUCTION 

The goal of this research is to develop theoretical, computational, and experimental tech- 
niques for predicting the effects of irregular topography on long range sound propagation in the 
atmosphere. Irregular topography here is understood to imply a ground surface that (l) is not 
idealizable as being perfectly flat or (2) that is not idealizable as having a constant specific acous- 
tic impedance. The interest of this study focuses on circumstances where the propagation is 
similar to what might be expected for noise from low-altitude air vehicles flying over suburban 
or rural terrain, such that rays from the source arrive at angles close to grazing incidence. 


PERSONNEL 

In addition to the. principal investigators, A. D. Pierce and G. L. Main, a graduate student, 
James Kearns, and two senior undergraduate students in mechanical engineering, Daniel Benator 
and James Parish, are presently working on the project. The students have up until now been pri- 
marily engaged in the construction of the experimental facility, in the construction of equipment, 
in the procurement of equipment and instrumentation, and in the testing of the components of 
the facility. Exploratory experiments are now beginning, with all three students participating. 
The theoretical work has up until now been carried out mostly by Pierce and Main, with tutorial 
sessions underway to develop Kearn’s participation in this phase of the research. 

All of the personnel concerned with the project visited NASA Langley Research Center in 
November 1985 and discussed complementary NASA and Georgia Tech research activities with 
the NASA technical officer, Dr. John Preisser, and his colleagues. 
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LABORATORY FACILITY AND INSTRUMENTATION 

The principal activity during the first year of the aubject grant haa been the conduction 
of a laboratory facility (Fig. 1) that will be uaed in subs^uent experiments. Major component 
of the laboratory are an acouatic aource, a model topographical surface with which the acoustic 
signal interacts, and a data acquisition and analysis system. The room housing the laboratory 
(Fig. 2) is a standard university small laboratory room of dimensions 6.2m by 6.1 m by 4 3 m 
The wall, are cinderblock , the Boor is tiled, and some walls are lined with shelves, so in no 
sense doe, this space approach the ideal of an anechoic chamber. However, the room, dubbed 
the -Atmospheric Sound Propagation Facility" within Georgia Tech, is solely dedicated to this 
project. The investigators have been developing the instrumentation system to be such that the 
echoes from walls, floor, and ceiling can be gated out in time. 

Base tables for topographic experiments 

Four tables were made to be used for the scale model experiments. Each table is 1.2 meter 
wide by 2.4 meter long and 0.9 meter high. The tables were constructed so that they could be 
bolted together, forming one table 4.9 meter long by 2.4 meter wide. The table top is CDX 
plywood, 2 cm thick. The table frame is constructed of two-by-six (5 cm by 15 cm) yellow pine 
grade #1 planks; the table legs are constructed of four-by.four (10 cm by 10 cm) yellow pine 

grade #2 beams The fasterners holding the table together are machine bolt, and wood screws, 
so the table assembly is fully portable. 

The table frame was made by running two 2.4 meter length twoby-sixes parallel to each 
other, 1.2 meter apart. Then five equally spaced two-by-sixe, were mounted in between these 
firs. two. Next, a shelf was cut into the four-b, -fours, so that they would fit into the corner, of 
the frame and still leave warn, of the frame resting on the shelf. Each four-by.four was bolted 

into ,h, frame with three machine bolts tor extra strength. Then the plywood was placed on the 
frame and secured with wood screws. 

In the experiments currently beginning, the four tables are bolted together and are being used 

m the 4.9m by 2.4 m configuration; the term •table’ implies this configuration the reminder 
of this report. 
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10 kV 

J, /.Spark 


Microphone 


Figure 1. General diagram of experimental configuration for studying sound propaga- 
tion over model topographical ridge. Here C is capacitor, R is resistor, P is 
power supply, IBM PC is IBM personal computer, ISC is RC Electronics A/D 
conversion instrumentation, A is amplifier, and P is preamplifier 








Figure 2. Photogr of interior of laboratory room used in the study, showing most of 
the relevan juipment and the experimental facility. 
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Laboratory scale model topographical ridge 

A curved surface (Fig. 3) was constructed to be mounted on the table and used as a 
laboratory scale model of a topographical ridge. The contour of the surface has the shape of an 
arc of a circle. 

Four basic templates (or ribs) for the ridge were cut from a piece cf 1.2 meter by 2.4 meter 
exterior plywood, 2 cm thick. The top edge of each template was an arc of a circle; the bottom 
edge was the chord of a circle. The chord was 2.4 meter long and the radius of the circle was 
such that the maximum height of the arc relative to the chord was 29 cm. Thus the radius of 
curvature of the arc was approximately 2.5 meter. 

Identical halves of the curved surface superstructure to the table assemblage were constructed 
as follows. For each half, a pair of templates were each secured to a 1.24 m by 2.4 m base board 
of CDX plywood (2 cm thickness) by nailing a strip of 3.8 cm by 3.8 cm yellow pine to each 
side of the template and then nailing the strips to the CDX plywood. The curved topographical 
surface was then achieved by bending plywood sheets, 0.5 cm thick and 1.24 cm wide, over the 
template arcs and then nailing the plywood to the arcs. 

Spark Generator 


A spark gap (Fig. 4 and Fig. 5) was constructed to serve as an impulsive acoustic source. 
Sound is generated when a sudden current surge occurs across a 2-3 mm air gap. As indicated 
in Fig. 6, a 10 kV power supply provides charge at the rated voltage to a 1 uf‘ capacitor. A 1 
Mf2 resistor is in series with the capacitor and the power supply; the voltage across the spark 
gap is virtually the same els that across the capacitor plates because of the negligible electrical 
resistance of the 0.5 cm diameter welding cables that carry current to and from the gap. The gap 
is between two 0.5 cm diameter tungsten electrodes that form the terminal points of the welding 
cables. The electrodes are held in position by a twopronged plexiglass fork which is mounted 
on a tripod stand (Fig. 7). A rotary grinder is used to shape and polish the electrode points. 
This spark gap generates an acoustical signal (Fig. 8) whose spectral content is dominated by 
frequencies of the order of 10 kHz and whose peak amplitude at a distance of the order of 1 m 
corresponds to roughly 110 dB. 
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Figure 3. Close-up photograph of oblique side view of curved ridge resting on table that 
was constructed for studying propagation effects of topographic.il ridges. 
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Figure 4, Design drawing of spark gap apparatus used in general, on of transient acoust.c 
pulses. 
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Figure 5. Photograph of spark jumping across gap between electrodes held in plexiglas 
frame. 
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spark gap 



Figure 6. D.agram of spark generation apparatus. Electrodes avai'able for use are either 
copper or tungsten. 
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Photograph of spark generation apparatus, showing tripod holding the poexi- 
■tlas frame with inserted electrodes. The power supply and capacitor are housed 
in a plexiglas box as a safety precaution. 
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Data Acquisition and Analysis System 

The data acquisition system is composed of microphones, amplifiers, an analog-to-digital 
converter, and an IBM personal computer (Fig. 9 and Fig. 10). The system is capable of gath- 
ering data at a rate of 500 kHz which can subsequently be processed by the PC. The acoustics 
laboratory VAX computer is in an adjoining room and is available for more extensive computa- 
tions and storage. The amplifiers were designed and built expressly for this project; other system 
components were purchased. In addition, an apparatus is being designed and constructed to 
quickly and precisely position the microphones at arbitrary points in the field. 

Bruel & Kjaer quarter-inch condenser microphones are used for making the necessary pre- 
cision sound pressure measurements. The microphone sensitivity (ratio of induced open circuit 
voltage to external acoustic pressure) for these microphones is certified by the manufacturer be 
be virtually constant for frequencies up to 70 kHz, so we expect them to yield a relatively undis- 
torted response to pulses predominantly composed of frequencies between 10 and 30 kHz. The 
microphones are linear in their response over a dynamic range of up to 180 dB with a sensitivity of 
0.1 mV/ubar. A Bruel ic Kjaer pre-amplifier and power supply are also part of the microphone 
assembly. The pre-amplifier has a very high inmput impedance and low parallel capacitance 
which are needed to maintain the flat frequency response. This high impedance is provided by 
a vacuum tube cathode follower at the input stage. The B k K 2801 power supply is used to 
provide voltage to the microphone and pre-amplifier. 

For the circumstances of the contemplated experiments the microphone assemblage open 
circuit voltage is typically of the order of 50 mV. Because such voltages are two low to exploit the 
full dynamic range of the analog-to-digital conversion instrumentation, two low current amplifiers 
were constructed to magnify the voltage signal. Each amplifier contains a Motorola LF351N 
FET operational amplifier microchip, which has a high voltage dew rate (13 V/us) and a flat 
response over a wide range of frequencies. The des.gn of the amplifier contains a non-inverting 
voltage amplifying circuit (Fig. 11). A variable gain is achieved by an array of feedback resistors 
controlled by an external multi-position switch. Frequency independent gains from unity up to 
100 are possible. The amplifier is powered by two parallel series of 9 volt batteries. The transient 
voltage level is guaranteed by parallel 0.1 uF capacitors. The amplifier is enclosed by a aluminum 
box which is grounded to the batteries. The box serves as a sh.eld against electromagnetic noise 
generated by the spark generator. Shielding considerations also motivated the use of a self- 
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contained power supply. 

The amplified analog signal is converted to digital form by an integrated hardware and 
software system produced by RC Electronics Inc. and called ”Computerscope ISC-16”. This 
system consists of a 16 channel A/D board which is inserted into a slot within an IBM PC, an 
external instrument interface, and the scope driver software. The system is capable of sampling 
data at a maximum aggregate rate of 0.5 MHz over as many as 16 channels. This state-of-the-art 
RC system is relatively new to the market and differs from other commercially available data 
interfaces for personal computers in its high data accession rate, which is adequate for acoustical 
experiments at frequencies in the 10’s of k.lohertzes range; the 0.5 MHz sampling rate provides 
50 data points per cycle for a 10 kHz signal. The system effectively transforms a personal 
computer into a low cost transient recorder or digital storage oscilloscope and should allow a 
greater flexibility in the digital processing of acoustical data. The system allows an input voltage 
signal with a peak to peak range of 20 volts centered at zero to be resolved to 12 bit accuracy 
or equivalently to 1 part in 4000. The incoming transient signal is stored within a 64 kilobye 
memory buffer. Various modes of triggering are possible. In particular, an external channel is 
provided exclusively for triggering without occupying any memory space, although it is possible 
to trigger off of another channel or off of threshold levels of slope or amplitude. The Scope Driver 
software allows for flexible manipulation and display of the captured data. The display is similar 
to that of an oscilloscope. We anticipate that all of the necessary spectrum analysis and transfer 
function calculations can be carried out, subsequent to data capture, by the host IBM PC, but 

the acoustic group’s VAX in the adjoining room is available for computations too involved oi 
lengthy for the personal computer. 
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Figure 9. Photograph of part of the experimental equipment showing microphone power 
supplies, the amplifiers whose design and construction are described in the text, 
and the monitor of the IBM personal computer. A typical transient acoustic 
pressure trace from a single microphone can be seen on the r<- •< <r screen. 




Figure 10. Photograph of the IBM personal computer with peripheral equipment which 
allows it to function as a digital storage oscilloscope or transient recorder. 
Corner of table facility of model topographical ridge can be seen at the far left. 
Cables from microphone assemblages lead to RC Electronics instrumentation 
interface, which in turn is connected to 16 channel high speed 12 bit A/D 
plug-in within the personal computer. 
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ANALYTICAL STUDIES 


Successful application of theoretical acoustics to outdoor propagation over undulating ter- 
rain is in principle possible, but presents challenges. The authors considerations are presently 
limited to when the terrain is slowly varying over distances comparable to a wavelength; many 
realistic outdoor situations should be well-modelled without violation of such a restriction. The 
overall hope is that asymptotic and matching techniques can enable one to splice together math- 
ematical models for intricate circumstances (such as multiple undulations) from those for simpler 
circumstances. 

Diffraction by a single ridge of finite impedance 

In the research program currently in progress, the understanding of diffraction by a single 
smooth ridge (Fig. 12) is a key element. Diffraction by a curved surface has a venerable and 
extensive literature, although much of it is specifically written for electromagnetic wave applica- 
tions. There is need for a readily assimilable treatment of acoustic diffraction by curved surfaces 
of finite impedance that is easily adaptable to servitude as a building block for a broader theory for 
propagation over irregular terrain. One desires simple analytical models or computational algo- 
rithms that are applicable on the surface and throughout the transition between illumination and 
shadow, not just deep within the shadow zone. Consequently, curved surface diffraction has been 
examined afresh, using the modern conceptual framework of matched asymptotic expansions. 

The theoretical work on the project to date has been especially influenced by the work of V. A. 
Fock [l|, who wrote a number of important papers on electromagnetic wave diffraction during the 
1940’s and 1950’3 that were later translated and republished together in a single volume. However, 
our method of derivation differs in some major details from that of Fock, and it is believed that 
the fresh perspective will facilitate the extension to broader classes of circumstances. 
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Creeping wave solution 


NAG-1-566 Semiannual Report, page 23 


Relatively simple results have been derived for the case when the source and listener are on 
opposite sides of a topographical ridge (Fig. 12) and the listener is deep into the shadow zone. 
With some miner distinctions, such results have been previously stated within an acoustical 
context by Hayek, Lawther, Kendig, and Simowitz [2], who adapt a theory developed by Keller 
[3] to diffraction by a cylinder-topped wedge of finite impedance. Our results, stated here for 
brevity without a derivation, extend those of Hayek et al. to cases where the radius of curvature 
and the surface impedance may vary with position. For simplicity, we consider source and listener 
to be on opposite sides rf the ridge; the extension to the oblique incidence case can be worked 
out without difficulty using relatively simple concepts [4]. 

The shortest path correcting source and listener has three segments, with lengths L lt L g , 
and L 2 . The segment of lei.gth L x is straight and terminates at the ridge at point a, where the 
segment is tangent to the curved surface. Similarly, the segment of length L 2 proceeds from a 
tangent point b to the listener position. The segment of length L g (< g for ‘ground’) proceeds along 
the curved top of the ridge from a to b. The surface’s local radius R(s) of curvature and specific 

impedance Z s (s) are functions of distance a along the surface. Here 3 = 0 corresponds to point 
a. 

The limiting case for which the simplest results most ideally apply is that where kL x , kL g , 
and kL 2 are all substantially larger than unity; here k is the wavenumber 2ir//c of the sound 
radiated by the source (strength 5). It is also implicitly assumed that the listener is well below 
the plane tangent to the ridge’s surface at point a, which separates the illuminated and shadowed 
regions. For this limiting case, the sound reaching the listener can be regarded as carried by a 
succession of ‘creeping waves’ that travel along the surface from point a and which shed rays into 
the shadow zone, each such ray proceeding along a straight line that if tangent to the surface; 
segment L 2 is a path of such a shedded ray. In the extreme limiting case of the sort considered 

above, the first creeping wave term dominates the sum and the complex amplitude of the acoustic 
pressure can be written 


P = 


Se ' 
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L'/* L \ /2 L l 2 /2 


R.i fib 
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Here L - L x + L 2 -h L g is total path length; /?„ and Rb are radii of 
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ground-induced phase shift 4> g and the ground-induced attenuation N a in nepers are respectively 


4 , = t/12+ f , (k/2R 2 ) 1 ^r R {q)ds 
J o 

( 2 ) 

", = [ L '(k/2R 2 y' 3 T l (q)ds 

J O 

( 3 ) 


Here t r and r, are real and imaginary parts of a quantity r that is the root having smallest 
imaginary part of the equation (prime denoting derivative) 


( 4 ) 


(<*) - qwi (a) = 0 

where tui(o) is a Fock function [1,5), given alternatively by 

= e’" /6 2* 1/2 Ai(ae°"/ 3 ) ( 5) 

in terms of the Airy function. The root a depends on a normalized surface admittance parameter 

q = i(kR/2)^pc/Z s (6) 

which varies with distance s if R or Z$ vary with 8. 

The remaining quantities A a and A b that enter into Eq. (1) are values at a and 6, respectively, 


of a quantity A(q ), defined such that 

A{q) = 2~ 5 / l2 jr~ 1 /* [ re -W 3 - 4 */®] ~ 1/3 [Ai(ae <2,r '' 3 )j _1 

T f the surface is rigid, then q = 0 and 

r = 1.0188e ,,r/3 = 0.5094 + i0.8823 

Ai(-re"’* /3 ) = 0.5357 
A{ 0) = 0.7817 

For small but nonzero q, an appropriate approximation is 

r = 1.0188e‘ ,t / 3 + a'"/ 0 g/1.0188 


( 7 ) 


(8a) 

(86) 

(8c) 


( 9 ) 


The correction affects the exponent factor N t and consequently may be of importance; in contrast, 
little harm is done if A{q) is approximated by A(0). A typical value of q can be estimated by 
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taking R — 3 m and / = 500 Hz. A survey of ground impedance data is given by Attenborough 
[6], who fits semi-empirical formulas to such data. Using his expression for the impedances of 
grassland reported by Embleton, Piercy, and Olson, one obtains 7.19 ■+■ «8.19 for Z s at 500 Hz. 
Such values then lead to an estimate of 


q = 0.22e*° 72 (10) 

Thus, one can regard q as small, but not necessarily negligibly small. 

Higher order terms in the creeping wave series have the same form as Eq. (1), the only 
distinction being that the calculations must use higher roots r„ of the transcendental equation 
(5), the roots being ordered according to the magnitude of their imaginary parts. The sum 
becomes slowly convergent, however, when L g becomes sufficiently small that many of the N g 
are close to zero. 
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Geometrical acoustics outer solution 
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At the top of the ridge and in the region of transition from illumination to shadow, a 
more nearly appropriate solution can be developed using the method of matched asymptotic 
expansions. A prototype twodimensional problem (Fig. 13) is when a plane wave of constant 
frequency with complex pressure amplitude P t exp (ikx) reflects and diffracts at a locally reacting 
(impedance Z s ) curved surface whose radius cf curvature R is not necessarily constant, but is 
nevertheless everywhere large compared with l/k. One argues with confidence that the field 
outside this surface for x < 0 can be satisfactorily predicted by geometrical acoustics [4]; this 
technique should also apply for sufficiently large positive y when x > 0. This general region 
is termed the outer region , because in the terminology- of matched asymptotic expansions, the 
geometrical acoustics solution for this region, when extrapolated down to the vicinity of the top 

of the surface (where y « 0 and x « R), furnishes the outer boundary condition for an inner 
solution that applies near the top of the barrier surface. 

The field in this outer region is a superposition of incident and reflected waves, such that 

P = P x e' k * + P % [A(Q) / A(l)\ l /*$ie ikx °e' ke ( U) 

where « is the reflection coefficient and A(£) denotes ray tube area after propagation a distance 

t from the reflection point. The reflection point (x 0 ,y 0 ), the local angle of incidence 9„ the local 

curvature radius R, a„u the reflected ray path length t can all be determined for given listener 

coo.dinates (x,y) using the law of mirrors and the the mathematical description of the surface 
(Fig- 14). 

Analysis of the so-derived geometrical acoustics solution for .he limiting case of points in the 
vicinity of the curved surfece’s top yields 


where 


P>e 


k 


’{■ 



1/2 




( 111 ) 


Q - [(4/9)x 2 4- (2/3) fly] 1 ' 2 
0 = (2 k/R 7 ){Q 3 - (8/27)x 3 - (2/3)rtxyj 


with R being the radius of curvature of the surface at the top (z = 0, y = 0). 
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Solution near top of ridge 

Scaling parameters L x and L y , equal to R/(kR ) 1/,a and R/^kRj 2 ^ 3 , can be introduced such 
that, when the pe ~' k * yielded by Eq. (12) above is expressed in terms of x/L x and y/L y , the 
resulting expression is independent of k and R. Since this furnishes the outer boundary condition 
on the inner solution, one anticipates that the inner solution should have comparable features. 

To develop the inner solution, the top of the surface is approximated by a parabola y = 
-x 3 /2 R, and the Helmholtz equation is expressed in parabolic cylinder coordinates u and v, such 
that 


X = u(l + [tl/ft]) 

(14a) 

y = v(l +[v/2R\)-u 2 /(2R) 

(146) 

so v = 0 corresponds to the diffracting surface. One then sets p equal to P t exp(tiku) times 
a function F of u/ L x and v/L y . The impedance boundary condition is Jso expressed in a 
nondimensional form using these variables. When the derivatives of F with respect to its nondi- 

mensionalized arguments are all regarded as being of the order of unity the terms 
differential equation satisfied by F become ordered by powers of (/c/Z)~ 1/3 . 

in the partial 

Substantial agreement with Fock’s notation is achieved if one set 3 


< = (2 /kR) l/ \ e = u/(2 l ' 3 L x ), n = 2 l/3 v/L y 

( 5 ) 

p= P ie ' k * e '< , / 3 G(£ >n>q ) 

(16) 

where 


q = t(kR/2) l/ * pc/Zs 

(17) 


is an appropriately scaled and nondimensionalized surface admittance. (Expected numerical 
values of q are discussed further below.) To lowest order in the expansion parameter e, the 
function G satisfies the parabolic equation 

idGjd$ + d 2 G/dr) 7 + r)G = 0 ( 18 ) 


with the boundary condition 


dG / dr i + qG = 0 at rj -- 0 


( 19 ) 
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The outer bounder, condition (her. imprecisely eta, ed, for brevity) ie the, E q . (16) match E q . (.2) 
at large positive y cr large negative x. 

The general eolation of the above posed boundary value problem can be developed by Fourier 
transform and complex variable techniques, with the result 

G (f,d„) = a- £[„«,-,)- 4^f^», ( o-n)]e-dn (20, 

where v(f) and «,,(() (as well as ,h, function, u( f ) and «.,( t ) defined further below) are Fock 
functions (S| and simply related to Air, functions of complex argument. (The precise definition 
of these functions is given further below.) The integral solution (20) is trivially related to what 
is termed (4,5) Fock', fc rm of the van dc, Po\-B,cmmc, diffraction formal,. 

Field on top surface of ridge 

One simple limiting case of interest is the acoustic prepare on the surface of the ridge, which 

1.9 


P = Ke' kx e' (,/3 C((,0,q) 

and corresponds to rj = 0. From Eq. (20) one obtains 


( 21 ) 


C(£,0, ? ) = x" 1 / 2 f°° r KjgMg) ~ ^(opwifa) . 

oo ( a ) “ <7^1 (a) 


;w -«-.£> (2i) 

However, the numerator in the bracketed term in the integrand here is a wronskian of two solutions 
O the A„y different,.! equation, so i, must be a constant. One finds, after plugging in the leading 
asymptotic expression for large positive r, ,ha, the constant „ simply 1, on, has 


V - w l t/ — 1 


and, consequently, 


(22) 


£(£A?) = ir ~ 1 2 J 


e % '* * 


da 


~ w [i a ) ~ 9^1 (a) ^ ( 23 ) 

Here „ is the normalised surface admittance defined in Eq. (17) and { i, Wllh 

u being interpreted a, being approximately ,h, distance , along the surface from the top of the 
r, g, down the -haded side. A more precise identification is tha, « is „ a „ lch 

the sum of the exponents „,e factor, s- and e'f- which appear Eq. (21). ,s ,ks. Thu, 
one would rewrite that equation as 


p= P,e’ kt G((, 0,q) 


(24) 
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A form of Eq. (23) that is more appropriate for computation at small to moderate values 
of £ ca-.i be developed by first deforming the integration contour from the real axis to a broken 
contour that goes from ooe 12 ’'/ 0 to the origin and then to oo along the positive real axis. If one 
uses 0 to denote distance from the origin along the first leg of this contour and recognizes that 


w 1 {0 e ' 7 ’ r '*) = e '*' 3 w 2 {0) (25 a) 

w\{0j 9 ’/ 9 ) = e’^V 3 0J) (256) 


one can derive 


G{Z,0, q) 



wi AP) ~ e' 2 "/ 3 qw 2 {0) d P + * 



„*a< 


u>i (a) - qw-(a) 


r da 


(26) 


The two integrals that appear here are highly convergent because at large positive real values of 
their arguments both and w 2 approach 


^(z) — u- 2 (z) 


1/-* c (2/3) 


(27) 


Thus one now has a version amenable to numerical computation. 

The apparent insertion loss 20 log(l/|G|) calculated using the above formula is plotted versus 
(; in Fig. 15. For the rigid barrier case, q = 0, the geometrical acoustics solution predicts a pressure 
doubling at the surface, so the insertion loss must approach -6dB at large negative £. 
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Transition between illumination and shadow 


A discussion of how the inner solution in Eqs. (16) and (20) above matches a further geomet- 
rical acoustics solution in the shadow zone is deferred to future reports. Deep in the shadow zone 
the appropriate version of the integral is a sum over residues from poles in the first quadrant, 

each such term giving rise to a creeping wave. This creeping wave solution is essentially the same 
as discussed earlier in the present report. 

The creeping wave series is not convergent at the boundary between illumination and shadow, 

where £ = r, 1 / 2 . An appropriate and suggestive form of the function G near this transition line 
when r/ i 3 somewhat larger than unuy is 






v'(s) - qe' ^ v(s) 


V w'(s) - <fe‘ V- u; 2 (a) 

U'(») -«».(■)/ 


da 


(28) 


where X (2/jt) tj I ({ r? 1/2 ) and A D (X) is the diffraction integral [4], which is simply related 
to Fresnel integrals and which is invariably present in asymptotic expressions for diffraction by 
sharp edges. Additional restrictions produce significant analytical simplification. 


Airy and Fock functions of complex argument 


As described in the preceding sections, the theory of diffraction by curved surfaces of finite 
acoustic impedance involves integrals (contour integrals in general) with integrands that can be 
expressed in terms of Airy and Fock functions of complex argument. It would therefore seem 
imperative that one have algorithms capable of calculating such functions to high precision for 
arbitrary complex argument. The algorithms we found reported in the literature were developed 
for functions of real argument only, so some efTort was devoted to developing new algorithms. 
The subroutines described in the present report are in a vers, on ( IBM Professional Fortran or. 
briefly. Profort) of FORTRAN 77 that can be used on the IBM PC, but it is intended that they 
be adapted in the near future to VAX Fortran. To achieve the desired accuracy with Profort 
(roughly ten significant figures), it was necessary to use double precision. IBM Professional 


if 


' i 


5 >. 




S?v'*” 
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Fortran (Profort) has complex number capabilities, but not in double precision, so the present 
version does all complex arithmetic explicitly. 

In principle, the Airy function of complex argument can always be calculated from the power 


series form 


where the leading coefficients are 


Ai(*) = cx/(z) - c 2 g(z) 


ci = 3" 2/3 /r(2/3) = 0.355028053887817 . . . 
c 2 = 3~ 1/3 /r(l/3) = 0.258819403792807 . . . 


and the intrinsic power series are 


^ 1 + 3-2* 6 • 5 • 3 • 2* + 9-8-6-5-3-2 z9 + "' 


g(z) = z + 


7 • 6 • 4 • 3 


10 • 9 • 7 • 6 • 4 • 3 


z lc> + ... 


In practice, however, this representation is useful for calculations only for moderately small 
arguments z. The program presented here uses it only if \z\ < 3. 

For larger arguments, one is initially tempted to use an asymptotic series representation lor 
the Airy function. However, such a series is not convergent absolutely. Although the magnitudes 
of successive terms may initially decrease, they eventually reach a minimum and then increase 
without limit. If one keeps only those terms up to and including the term of minimum magnitude, 
then this is as good as one can do with an asymptotic series. The error is of the order of magnitude 
of the next neglected term. Some trial calculations suggested this would not be good enough 
(given a desired precision of at least 1 part in 10°) when !*| was of the order of 3 Consequently, 
an alternate procedure was used. 

To describe this alternate procedure, one first notes that the Airy function can alternately 
be described by the contour integral 


.« I • / 3 + t » 


The contour C Kx can be initially thought of as proceding in the complex s-plane along the broken 
line which goes from ooe ,5,r/a to the origin and then to oce”' / '’. If the argument variable z 
lies, however, in the right 2/3-rd’s of the complex z-plane, then the integration path C A , can 
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be deformed to one that passes through a saddle point, going to this saddle point up a path of 
steepest descents and then away from this saddle point down a path of steepest descents. (By 
the statement that the argument variable z lies in the right 2/3-rd’s of the complex plane, one 
means that the phase of * lies between -2 jt/ 3 and 2*/3. Because of the identity in Eq. (47a), 
given further below, this turns out to be no real restriction.) 

The applicable saddle point, obtained by setting the derivative of the exponent to zero, is at 
9 = xz t ■ To change the contour to the path of steepest descents, one sets s = tz 1 ^ 2 +u such that 
the exponent in the integrand can be written -(2/3 jz 3 '' 2 - O where f? = z^u 2 - (»/3)u 3 . The 
saddle point now corresponds to u = 0 or, equivalently, to t = 0. The path of cteepest descents 
is a path along which l is real; the definition of i can be refined such that the mapping of C Ai to 
the 4-plane can be deformed to a path that goes from -oo to oo along the real axis. The integral 
expression for the Airy function can accordingly be rewritten 


Ai(z) = ^-e < 2 /3>* ( 3/2; J 
where t and u are related by 


21 

2UZ 1 / 2 - iu 2 


.-r 


dt 


e = z^'v - (i/3)u 3 


(33) 


(34) 


The latter is a cubic equation for u as a function of t\ the desired root must be zero when 
t is zero; moreover, u(t) must be a continuous function of i. The two possibilities correspond 
to u{t) for small l being either +l/z l/ * or -t/z l/i . The requirement that a contour from -oo 
to oo along the real t axis be an admissible deformation of the mapping of C Al into the f-plane 
indicates that the former choice is correct. Thus one can write 


u = 


Kl 

e 1 / 4 


(35) 


where h [l) is such that A( 0) — 1, and is a solution of the cubic equation 

i = *>-•«£ 

3 z 3 /< 

The appropriate solution of the above cubic equation can be worked out with some effort, 
the result being 


xz 


3/4 


* • 


where 


1 + e'* n A f e~' n/3 A- 1 } 

(37) 

ZP si;. 3‘^1 2 /3 

4 z 3 / 2 2 z 3 / 4 . 

(38) 


i 


V 


h 

i 

i 
i 
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For small values of lz~ 3 '* it is appropriate to replace the latter by its power series, which to 
fourth order is 

A-l +jj. +... (39) 

where here we abbreviate s = (3 l/2 /2 )iz~ 3 !*. 

the integral of e over t from -oo to oo is a -1 / 2 , one can derive from the above 
expressions 

e -(2/3), J '> r 

A • / \ C 


where 


, 3 /«\ - 


j F M {l,z 3 ^)e~ t7 dt j 

(40) 

• lf 3 ^-i , 

* 2z 3 / 4 1 

(41) 


or, equivalently, with the substitution of Eq. (37), 


F (P , 3 /<> - ~ 2 - 1 - e ia "/*A 2 - e -* 2 /p</3 /l -2 
e* 2,, / 3 A 3 + e~' 2 "/ 3 A~ 2 -)- 1 


with 


A ±2 = f(l + ^L) 1 /' 2 T ^ 

l V 4-3/2/ T 0,3/4 


3 l/2£ 


4/3 


(42) 


(43) 


4-2 3 / 2 ' ' 2z 3 / 

The leading term in Eq. (40) (i.e., that which results when F M is formally set to zero) is the first 
term in the asymptotic expression for Ai(z). 

What is achieved with the introduction of Eq. (40) is that the integrand is not oscillatory, 

so the integral is highly convergent. The integral is don. numerically us.ng a Hermite integration 
scheme [7], so that 


/ OO 10 

Fm * 3/4 ) U = Y. ** e ~ e ’ l f M(l,* 3/4 ) + F M (-ei, *3/4)] 

4=1 


(44) 


The sampling points t, and weights w, are tabulated in the listing of the subroutine aamairy. 

A similar procedure has been derived for computation of the derivative of the Airy function, 


AiTzl = --i!H,-i2/3). ,/ ’ f. | 1 r r it 3 / 4 \ ~ 

V ' 2x‘/ J l 1 ^TTI / G «(^ ' )« 

J — OO 


where 


G u (t,z 3 '*) = _i _ 2> < g '* 73 A 4 - e -*"/ 3 A ~ 1 

z 3/4 j + Jln/Stf + e -t2n/3 A -J 


dt 


(45) 


( 46 ) 


* 
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As mentioned previously, these above integral expressions, Eqs. (40) and (46), are valid 
only if the phase of z lies between -2 jt/ 3 and 2 jt/ 3. However, one can use these expressions in 
conjunction with the Airy function identities 

Ai(z) = e^AiCze-* 2 */ 3 ) + e - tw ' 9 Ai(xe- i49 ' 3 ) (47a) 

Ai'(z) - e _ ’’ r / 3 Ai'(ze -<2 ’ r / 3 ) + e**/ 3 Ai'(ze- 4 "/ 3 ) (476) 

Note that, if the phase of z is between 2z/3 and 4ir/3, then the arguments ze - ’ 2 */ 3 and ze— i4x/3, 
which appear on in the terms on the right sides of the above two equations, have phases between 
-2x/3 and 2x/3; thus each such term can be calculated using Eqs. (40) and (46). 

For the computation of the Fock functions and their derivatives, we use the relations 



v(z) = * 1/2 Ai (z) 

(48a) 


«>i(z) = e ,,r/a 2jr 1/2 Ai(ze* 2 " /3 ) 

oo 


t* 2 (z) = e-’ > / 6 2x 1 / 3 Ai(ze- 2 ^ 3 ) 

(48c) 

such that 

t/(z) = T 1/2 Ai'(z) 

(49a) 


w[(z) =e tSn / 6 2z 1 / 2 Ai , (ze >2 ’ r / 3 ) 

(496) 


«4(z) = e-^^x^Ai'tze-' 2 "'' 3 ) 

(49c) 

The core algorithms are 
argument. 

consequently those that evaluate Ai(z) and Ai ; (z) for arbitrary complex 

The algorithms given here have been checked against Fock’s tables (which appear on pages 
393-412 of his Electromagnetic Propagation and Diffraction Problems [l]). Fock tabulates u(z). 

u'(z), v(z), and v'(z) for real z between -9 and +9 to four significant figures. The function t>(z) 
is just x I/2 Ai(z), while 


u(z) = i(wifz) -f u/ 2 (z)) 

(50) 

such that 

Wi(z) = u(z) + iv(z ) 

(51a) 


w 7 (z) = u(z) - iv(z) 

(516) 
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If z is real, then both u(z) and t/(z) are real, u(z) and v(z) are the real and imaginary parts of 
u>i(z), or, equivalently, the real and negative imaginary parts of tu 2 (z). Our program’s results 
(believed to be accurate to 10 significant figures), when rounded off to 4 figures, agree identically 
with Fock’s results. For example, Fock gives u' (9) = 113. lOx 10 6 , and we find it to be 113. 095831 x 
10 *. 


PAPERS AND PUBLICATIONS 

The following paper will be presented at the forthcoming meeting of the Acoustical Society 
of America in Cleveland, Ohio in May 1986. 

Curved surface diffraction theory derived and extended using the method of 
matched asymptotic expansions. Allan D. Pierce, Geoffrey L. Main, and James A. 
Kearns, School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, Geor- 
gia 30332. — Consideration is given to the top of a wide barrier with variable radius of 

curvature R . The surface haw finite acoustic surface impedance Z . Because kR is assumed 
large, the illuminated region can be approximated by geometric acoustics, such that plane 
wave reflection rules apply locally. The intricate interference pattern between incident and 
reflected ray fields assumes a tractable analytical form near the barrier top, which is sub- 
sequently used in a MAE solution of the overall diffraction problem. Unambiguous length 
scales result for radial and tangential distances along the barrier top. The inner solution is 
developed by expressing the wave equation in terms of such scales, subsequently identify- 
ing the expansion parameter as (A:i?) -1 / 3 . A parabolic equation emerges, with a boundary 
condition involving a scaled impedance (Z/pc)(kR)~ l ^\ the outer boundary condition re- 
sults from matching to the geometric acoustics solution. Outer expansion of the solution of 
the parabolic equation into the shadow zone yields an inner boundary condition on the ray 
theory solution for the diffracted wave. Results are similar to those previously derived for 
electromagnetic diffraction problems by V. A. Fock, but the MAE interpretation facilitates 
an extension to problems of multiple barriers. (Work supported by NASA-Langley Research 
Center.) 
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The following paper will be presented at the forthcoming International Congress of Acoustics 
in Toronto in July 1986, and will appear in the proceedings of that congress. 

Sound propagation over large smooth ridges in ground topography. Allan D. 
Pierce, Geoffrey L. Main, James A. Kearns, Daniel R. Benator, and James R. Parish, Jr. 
School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332. 

A theory similar to those developed by Fock and others during the 1940’s and 1950’s for 
electromagnetic wave diffraction by curved surfaces applies to acoustic propagation at low 
angles with the ground over an intervening ridge of finite impedance. The creeping wave 
series is not used at the top of the ridge or for the transition between illumination and 
shadow; the analysis reduces instead to numerical and approximate integration of Fock’s 
form of the van der Pol-Bremmer diffraction formula. Laboratory scale experiments are in 
progress to test and guide the analytical developments. 

The following paper will be presented at the forthcoming 1986 International Conference on 
Noise Control Engineering (Inter-Noise 86) in Cambridge, Massachusetts in July 1986. 

Sound propagation over curved barriers. Allan D. Pierce, Geoffrey L. Main, James A. 
Kearns, and H.-A. Hsieh, School of Mechanical Engineering, Georgia Institute of Technology, 
Atlanta, Georgia 30332. — A general discussion is given of wide barriers with curved tops; 
examples of such are naturally occurring topographical ridges and earth berms with rounded 
tops. The analytical developments reviewed are for circumstances when the local radius of 
curvature R of the barrier is continuous along the surface and large compared to a wavelength. 
If the source and listener are at large distances from the barrier top and the listener is deep 
within the shadow zone, then the creeping wave series previously introduced into noise control 
applications by Hayek and others gives simple and accurate predictions. The present paper 
extends this model to instances where the surface impedance varies with position along the 
surface. The latter part of the paper introduces a matched asymptotic expansion theory 
that contains concepts and results analogous to those developed by V. A. Fock. Explicit 
numerical results are given for the acoustic pressure on the surface of the barrier near the 
point where acoustic shadowing begins. The extended theory also yields simple- results for 
the farfield transition between illumination and shadow. 
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APPENDIX — COMPUTER PROGRAMS 


Listing of computer programs for Airy and Fock functions 


The input and output subroutine, given here are temporary and intended only tor checking 
out the algorithm, with a desk-top monitor. The actual program per se consists of subroutine Airy 
and all those subroutine, that it call.. Subrouti- Fork uses Airy to compute the Fock functions. 
For the numerical evaluation of integral, that describe the diffraction by curved surfaces of finite 
impedance, programs will be written that call these two subroutines. 


program Airychek 

* Allan D. Pierce 

12/25/85 

double precision x,y,airyr,airyi,dairyr,dairyi, 

+ vr ,vi,dvr,dvi,wlr,wli,dwlr,dwli, 

+ w2r,w2i,dw2r,dw2i,r, pi, angle 

call input ( t, angle) 

pi = 3.1415926535897932D0 

x = r*dco«(angle*pi/130DO) 

y = r*istn(angle*pi/180D0) 

call 4irj/(x,y,airyr,airyi ) dairyr,dairyi) 

call Foci(x,y,vr,vi,dvr,dvi,wlr,wli,dwlr,dwli, 

+ w2r,w2i,dw2r,dw2i) 

call prtntAtr (x,y,airyr,airyi,dairyr,dairyi, 

+ vr,vi,dvr,dvi,wlr,wli,dwlr,dwli, 

+ w2r,w2i,dw2r,dw2i) 

step 

end 


subroutine mput(r, angle) 

double precision r, angle 
write (*,*) ’Program Airychek’ 
write (*,*) ’Version of December 1985’ 
write (*,*) 
write (*,*) 

write (*,*) ’What is magnitude of argument? ’ 
read (*,*) r 

write (*,*) ’What is phase angle in degrees? ’ 
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read (*,*) angle 

return 

end 


subroutine Foci(x,y 1 vr,vi,dvr,dvi,wlr,wli,dwlr,dwli, 
+ w2r,w2i,dw2r,dw2i) 

double precision x,y,airyr,airyi,dairyr,dairyi, 

+■ vr,vi,dvr,dvi,wlr,wli,dwlr,dwli ( 

+ w2r,w2i,dw2r,dw2i,pi,ar,ai,cr,ci , 

+ er,ei ) xl ) yl,x2,y2,sr,si,tr ) ti 

pi = 3.1415926535897932DO 
ar = daqrt( pi) 
ai = 0D0 

call Airy(x,y,airyr,aiiyi,dairyr,dairyi) 
call cprcx^ar.ai.airyr.airyi.v^vi) 
call eprod(ar,ai,dairyr,dairyi,dvr,dvi) 
cr = d:oa[ pi/1.5D0) 
ci = datn(pi/1.5D0) 
call cprod(x,y,cr ) cipcl,yl) 
er = dcos(pi/6D0) 
ei = dsin (pi/6D0) 
ar = 2D0*ar 

call 4»ry(xl,yl,airyrairyi,dairyr,dairyi) 

call cproa^ar.a^er^i^r^i) 

call cpro<f(sr,si,airyr,airyi,wlr,wli) 

call cpro<f (sr,3i,cr,ci,tr ,ti) 

call cprod(tr ) ti,dairyr,dairyi ) dwlr,dwli) 

call cprod(r, y,cr,-ci,x2,y2) 

call /liry (x2,y 2 ,airyr,airyi, dairy r, dairy i) 

call cprod (sr, -si, airy r, airy i,w2r,w2'.) 

call cprod(tr,-ti,dairyr,dairyi,dw2r,dw2i) 

return 

end 


subroutine A tr y(x ,y,airyr, airy i, dairy r, dairy!) 

double precision x,y,airyr,airyi,dairyr,dairyi, 
+ pi,cl,c2/ork,r 

integer N 

data pi /3.1415926535897932DO/, 

+ cl / 0.3550280538878 17 D0/, 
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+ c2 / 0.258819403792807 DO/, 

+ fork / 3.0DO/, 

+ N / 20/ 

r = dsqrt(x**2-|-y**2) 

If (r .le. fork) 

+ then 

call i4»ryl (x,y,airyr,airyi ) dairyr,dairyi, 

-t- N,cl,c2) 

else if (r .ge. fork) 

+ then 

call Airy2 (x,y, airyr, airyi, dairyr, dairyi, pi) 

end if 

return 

end 


subroutine Airy 1 (x,y,ahyr, airyi, dairyr, dairyi, 
+ N,cl,c2) 

double precision x,y,airyr,airyi,dairyr,dairyi, 
+ br,bi,cr,ci,zetr,zeti,fr,fi, 

+ gr 1 gi,dfr ) dfi,dgr ) dgi,cl,c2 

integer N 


br = x 
bi — y 

call cprod{x, y,br,bi,cr,ci) 
call cprod{x, y,cr,ci,zetr,zeti) 
call 3eratry(zetr,zeti,N, frfi.gr, gi, 
+ dfr,dfi,dgr,dgi) 

call cprod (x,y,gr,gi,br,bi) 
airyr - cl*fr - c2*br 
airyi = cl - c2*bi 
call cprod(cr,ci,dfr,dfi,br,bi) 
dairyr = 0.5D0*cl*br - c2*dgr 
dairy i = 0.5DO*cl*bi - c2*dgi 

return 

end 


subroutine Airy 2 (x,y, airyr, airyi, dairyr, dairy i ,pi ) 
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double precision x,y,airyr,airyi,dairyr,dairyi, 

+ pi, r,phi,ar,ai,incr,inci, phase, 

+ ara,aia,dra,dia,arb,aib,drb,dib, 

+ phia,phib,phic,arc,aic,drc,dic, 

+ uar,uai,ubr,ubi,ur,ui,ci 

r = dsqrt(x** 2 + y**2) 
phi = pAase(x,y,pi) 

if (phi .gt. (4D0/3D0)*pi - ID-14) phi = phi - 2D0*pi 
if (phi .It. pi/1.5D0 - ID-14 .and. 

+ phi .gt. - pi/1.5D0 + ID-14) 

+ then 

call Ai'ry2a (r,phi,airyr,airyi,dairyr, 

+ dairyi.pi) 

else if (phi .gt. pi/1.5D0 + ID-14 .or. 

+ phi ,!t. -pi/1.5D0- ID-14) 

+ then 

phia = phi - (phi/dabs(phi))*pi/1.5D0 
call Airy 2a (r,phia,ara,aia,dra,dia,pi) 
phib = phia - (phi/dabs(phi))*pi/1.5D0 
call Airy2& (r,phib,arb,aib,drb,dib,pi) 
ar = 0.5D0 

ai = (phi/ dabs(phij) *dsqrt(3D0) *0.5D0 

call cprod ( ar,ai,ara,rda, airyr, airy i) 

call cprod (ar,-ai,dra,dia, dairy r, dairy i) 

call cprod(ar,-ai,arb,aib,incr,inci) 

airyr = airyr + incr 

airyi = airyi 4- inci 

call cprod (ar ,ai,drb, dib, incr , inci) 

dairyr = dairyr + incr 

dairyi = dairyi -f inci 

else 

phia = pi/1.5D0 - ID-15 

call Airy 2a (r,phia,ara,aia,dra,dia,pi) 

phib = - pi/1.5D0 + ID-15 

call A\ry2& (r,phib,arb,aib,drb,dib,pi) 

phic = 0D0 

call Airy2a. (r,phic,arc,aic,drc,dic,pi) 
ar = d«?rt(3DO)*0.5DO 
ai = 0.5D0 

call cprod (ar,ai,ara,aia,uar,uai) 

call cprod (ar,-ai,arb,aib,ubr,ubi) 

ur - (uar -f ubr)/2D0 

ui = (phi/ da6»(phi))*arc/2DO 

ci = (phi/ da6s(phi))*ai 

call '•prod(ar,-ci,ur,ui,airyr,airyi) 

call cp rod (ar ,- ai , dra, dia, uar ,uai) 


•o-.nr- « 
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call cpr£>d(ar,ai,drb,dib,ubr,ubi) 
ur = (uar + ubr)/2D0 
ui = - (phi/ <fa6a(phi))*drc/2DO 
call cprod (ar,ci,i>r ,ui,dairyr,dairyi) 


end 

return 

end 


subroutine Airy2a (r,phi,airyr,airyi ) dairyr,dairyi,pi) 

double precision r ) phi,airyr,airyi ) dairyr,dairyi,pi, 
+ U,el,e2,a,ar,ai,cr,ci, 

+ b,br,bi,fr,fi,dfr,dfi 


u = r*dsqrt{r) 

el = u*</cos(l.5*phi)/1.5DO 

e2 = u*d«'n(1.5*phi)/1.5D0 

a = dexp(.e 1) / (2D0*<fo ? rt(pi*d« ? rt(r))) 

ar -- &*dcos(e2 +0.25D0*phi) 

ai = - a*<fsi«(e2 +0.25D0*phi) 

call aamatry(r,pki,fr,fi,dfr,dfi,pi) 

call tprod(ar,ai,fr,fi,cr,ci) 

airyr = cr 

airyi = ci 

b = d3qrt(daqrt(T) /pi)* dcxp(-el)/2D0 
br = - b*deos(0.25D0*phi - e2) 
bi = - b*dsm(0.25D0*phi - e2) 
call cprod(br,bi,dfr,dfi,cr,ci) 
dairyr = cr 
dairyi — ci 


return 

end 
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function pAa$e(x,y,pi) 

double precision x.y.px.py.phase.r.pi 
if (x eq. 0D0 and. y eq. 0D0) 
t- then 

phase— 0D0 
return 
end if 

if (fa6s(x) ge. dabs(y)) * 
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4- then 

py = das in( dabs(y)/ dsqrt(x**2+y**2) ) 
if (x .ge. 0D0 .and. y .ge. 0D0) 

4- then 

phase = py 

else if (x .le. ODO .and. y .ge. 0D0) 

4 - then 

phase = pi - py 

else if (x .le. ODO .and. y .le. ODO) 

4 - then 

phase = pi + py 

else if (x .ge. ODO .and. y .le. ODO) 

4- then 

phase = - py 
end if 

else if ( dabs(x ) .le. dabs(y)) 

4 - then 

px = dasin( dabs(x)/ dsqrt(x** 2+y**2) ) 
if (y .ge. ODO .and. x .ge. ODO) 

4- then 

phase = 0.5D0*pi - px 
else if (y ge. ODO .and. x .le. ODO) 

4 - then 

phase = O.5D0*pi 4- px 
else if (y .le. ODO .and. x le. ODO) 

+ then 

phase = 1.5D0*pi - px 
if (phase .gt. 4DO*pi/3DO) 

-I- then 

phase =; phase - 2DO*pi 
end if 

else if (y .le. ODO .and. x .ge. ODO) 

4 then 

phase = -0.5DO*pi px 
end if 
end If 

return 

end 



4 


€ 


subroutine cprod (ar,ai,br,bi,cr f ci) 
double precision ar,ai,br,bi,cr,ci 


cr — ar*br - ai*bi 

ci ~ ai*br + ar*bi < 






f 
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return 

end 


subroutine seratry (zetr ,zeti,N,fr ,fi,gr ,gi ,dfr ,dfi , 
+ <*gr,dgi) 


double precision zetr.zeti,fr,fi,gr,gi dfr.du, 
+ dgr,dgi,denom 

integer Nj 


fr = IDO 
fi = 0D0 
gr = IDO 
gi = 0D0 
dfr = IDO 
dfi = ODO 
dgr = IDO 
dgi = ODO 


do 30 k - 1,N 
‘ j = N +1 - k 


denom = (3DO*j)*(3DO*j-lDO) 
call oneatp(denom,zetr,zeti,fr,fi) 
denom = (3D0*j+lD0)*(3D0*j) 
call onestp{denom,zetr,zeti,gr,gi) 
denom = (3D0*j+2D0)*(3D0*j) 
call or»eaip(denom,zetr,zcti,dfr,dfi) 
denom = (3DO*j)*(3DO*j-2DO) 
call onestp (denom, zetr ,zeti, dgr, dgi) 
30 continue 

return 

end 


subroutine o-\estp (denom, zetr, zeti.fr ,fi) 

double precision denom, zetr, zeti,fr,fi, 

+ br,bi,cr,ci 


br = zetr/denorn 
bi = zeti/denom 
call cprod(br,bi,fr,fi,cr,ci) 


i 






NAG-1-566 Semiannual Report, page 48 



I 

ii 


i 

! 

I 

r 

i 

i 

n 


subroutine asmairy (r,phi,fr,fi,dfr,dfi,pi) 


double precision z(10), w(10) 

double precision intr, inti, dintr, dinti.fr ,fi,dfr, 

+ dfi,r ,phi,phia,eta,phieta,u,pi, 

+ krl,kil,dkrl,dkil,kr2,ki2,dkr2,dki2 

data z(l) / 0.2453407083009 DO/, 

+ z(2) / 0.7374737285454D0/, 

+ z(3) / 1.2340762153953DO/, 

+ z(4) / 1.7385377121 166 DO/, 

+ z(5) / 2.2549740020893D0/, 

+ z(6) / 2.7888060584281D0/, 

+ 2(7) / 3.3 178545673832D0/, 

+ 2(8) / 3.9447640401 156D0/, 

+ 2(9) / 4.6036824495507D0/, 

+ 2(10) / 5.38748089001 12D0/ 

data w(l) / 0.49092 15006667DO/, 

+ w(2) / 0.493843385272 1D0/, 

+ w(3) / 0.49992087 13363 D0/, 

+ w(4) / 0.5096790271 175DO/, 

+ w(5) / 0.5240803509486DO/, 

+ w(6) / 0.54485 17423644DO/, 

+ w(7) / 0.5752624428525DO/, 

+ w(8) / 0. 622278696191 4D0/, 

+ w(9) / 0.704332961 1769D0/, 

+ w(10) / 0. 89859196 14532 D0/ 


intr = 0D0 
inti = 0D0 
dintr = 0D0 
dinti — 0D0 


phia = phi 

if (phia ,gt. pi) phia - phia-2.0D0*pi 
if (phia .It. -pi) phia = phia+2.0DO*pi 
eta = i«?rt(dsqrt(r*r*r)) 
phieta = (3D0/4D0)*phia 


do 20 j — 1 , 10 
u = z( j)/eta 
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t 
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call steep (u.phieta.pi.krl.kil.dkrl ,dkil ) 

call «tee/>(-u,phieta,pi,kr 2 ,ki 2 ,dkr 2 ,dki 2 ) 
mtr = intr + w(j)*(krl+kr2)*dezp(-z(j)*z(j)) 

inti = inti + w(j)*(kil+ki2)*rfexp(-z(j)* z (j)) 
dintr = dintr + w(j)*(dkrl+dkr2)*dexp(-z(j)* z (j)) 
dinti = dinti + w(j)*(dkiH-dki2)*dex P (.z(j)* z (jl) 
20 continue 

fr = IDO + intr/ dagrt(pi) 
fi = inti/ dsqrt (pi) 
dfr = IDO •+■ dintr/ dsqrt( pi) 
dfi = dinti/ <fo^r£(pi) 

reixrn 

ena 


subroutine at«p(u,phieta,pi,kr,ki,dkr,dkj) 

double precision u,phieta,pi,kr,ki ( dkr,dki, 
ub.ubr.ubi.dumr.dumi.ubsqr, 
ub3qi,ucubr,ucubi,rr,ri,denom, 
radsqr .radsqi /ad ,pt , phase , 
phirad,a32r ) a32i,a32,pa32, 
asq, phasq,a,pha, denomr ,denomi, 
denomsq,recipr,recipi,yr,yi, 
numr^umi.newr^ewijXrjXi^r^i 

kj 

ub = *?rt(3D0)*0.5DO*u 
ubr = ub*dcoa(phieta) 
ubi = ■ ub*dsm(phieta) 
dumr = ubr 
dumi = ubi 

call cprod(dumr , dumi, ubr, ubi, ubsqr.ubsqi) 

if (da6s(ub) .It. O.OOIDO) 

■+■ then 

call cprod(ubr,ubi,ubsqr,ubsqi,ucubr.ucubi) 
numr = (8D0/9D0)*ubr - (40D0/243D0)*ucubr 
numi = (8DO/9DO) *ubi - (40D0/243D0) *ucub. 
numr = numr - * ? rt(3D0)*(14D0/8!D0)*ubsqi 
numi = numi + <h ? rt(3DO)*(14DO/81DO)*ubsqr 

denomr = - numr 

denomi = - 4D0/ dsqrt{ 3D0) - numi 
recipr = denomr / (denomr**2 -f denomi**2) 


+ 

+ 

+ 

+ 

+ 

+ 

+ 
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i 
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recipi = -denomi / (denomr**2 + denomi**2) 
call cprod(recipr,recipi,numr,numi,kr,ki) 
newr = (32D0/9D0)*ubr - (100DO/243DO)*ucubr 
newi = (32D0/9D0)*ubi - (10ODO/243D0)*ucubi 
newr = newr + dsqrt(Z DO) * ( 10DO/ 8 1 DO) *ubsqi 
newi = newi - d*grt(3DO)*(10DO/81DO)*ubsqr 
call cproc^recipr^ecipijnewrjiewi^dkrjdki) 
else if (dabs(ub) .ge. O.OOIDO) 

-t then 


i 


i 


j 


\ 


t 

i 


\ 

i 


; 


| 




if (ub .gt. 10D0) 

+ then 

rr = ubsqr/(ubsqr**2+ubsqi**2) 
ri = -ubsqi/(ubsqr**2-|-ubsqi**2) 
a32r = 0.5DO*ubr/(ubr**2 + ubi**2) 
a32i = -0.5D0*ubi/(ubr**2 + ubi**2) 
do 20 j = 1,10 
k = 5 - j + 1 

denom = (0.5D0 - k)/(k + 1D0) 
call onestp(denom,rr,ri,a32r,a32i) 

20 continue 

else if (ub .le. 10D0) 

4 - then 

radsqr = IDO 4- ubsqr 
radsqi = ubsqi 

rad = dsqrt(dsqrt( radsqr**2 4- radsqi**2)) 
pt = p/»ase(radsqr, radsqi, pi) 
if (pt gt. pi) pt = pt - 2D0*pi 
if (pt It. -pi) pt = pt 4- 2D0*pi 
phirad = pt/2DO 
a32r = rad*dcos(phirad) - ubr 
a32i = rad *dain (phirad) - ubi 
end if 

a32 = dsqrt{ a32r**2 4- a32i**2) 

pa32 = phase( a32r,a32i,pi) 

if (pa32 gt. pi) pa32 = pa32 - 2D0*pi 

if (pa32 .It. -pi) pa32 = pa32 4- 2D0*pi 

asq = a32**(4DO/3DO) 

phasq = pa32*(4DO/3DO) 

a = di^rt(asq) 

pha = phaaq/2D0 

denomr = IDO 4 - (asq+lDO/asq)*dcos(phasq-<-pi/1.5DO) 
denomi — (a9q-lDO/asq)*d«m(phasq+pi/l .5 DO) 
denomsq = denomr**2 + denomi**2 
recipr ~ denomr/denomsq 
recipi = -denomi/denomsq 
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yr = -2D0*u*d«t'n(phieta) 
yi = -2D0*u* dcoa (phieta) 
numr = yr - denomr 
numi = yi - denomi 

call cprod(recipr,recipi,numr,numi,kr,ki) 
xr = (a + 1 DO/a) * dcoa (pi /3D0 + pha) 
xi = (a - lD0/a)*<fo*n(pi/3D0 + pha) 
call cprod(yr,yi,xr,xi,br,bi) 
numr = br - denomr 
numi = bi - denomi 

call cprod(recipr,recipi,numr,numi,dkr,dki) 
end if 

return 

end 
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subroutine print Air (x,y,airyr,airyi,dairyr,dairyi ) 
+ vr,vi,dvr,dvi,wlr,wli,dwlr,dwli, 

+ w2r,w2i,dw2r,dw2i) 

real*8 x,y, airyr, airyi, dairyr,dairyi,vr,vi, dvr, dvi, 

+ wlr,wli,dw Ir,dwli,w2r,w2i,dw2r,dw2i 

write (*,*) 
write (*,*) 

write (*,1) ’ x = x, ’y = y 
1 format (10X, A, F10.4, 10X, A, F10.4) 

write (*,*) 

write (*,3) ’Airyr’, ’Airyi’, ’Derivr’, ’Derivi’ 

3 format (12X,A,11X,A,10X,A,10X,A) 

write (*,*) 

write (*,101) airyr, airyi, dairyr, dairyi 
write (*,*) 

write (*,3) ’ vr’, ’ vi’, ’ dvr’, ’ dvi’ 
write (*,*) 

write (*,101) vr,vi,dvr,dvi 
write (*,*) 

write (*,3) 1 wlr’, wli’, ’ dwlr’, ’ dwli’ 
write (*,*) 

write (*,101) wlr, wli, dwlr, dwli 
write (*,*) 

write (*,3) ’ w2r’, ’ w2i\ ’ dw2r’, ’ dw2i' 

write (*,*) 

write (*,101) w2r,w2i,dw2r,dw2i 
101 format (lx, 4D16.8) 





£ 


i 


i 


* 


G 



return 


end 
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Program to compute pressure on ridge surface 


Here we list the computer program that was used in the computation of the curves shown 
in Fig. 15 for the apparent insertion loss of the topographical ridge, 
program ridgchek 


Allan D. Pierce 
2/13/86 

double precision xi,qr,qi,pi,gr,gi,g,lo8s 

call ridgput (xi,qr,qi) 
pi = 3.1415926535897932D0 
call n'dytn^xi.qr.qi.gr.gi.pi) 
g = daqrt{ gr**2 4- gi**2) 
loss = 20D0* dlog 10( IDO /g) 
call prtridge(xi, qi ,qi ,gr ,gi g,loss) 

stop 

end 


subroutine ru/ 0 tn£(xi,qr 5 qi,gr,gi,pi) 

precision xi,qr,qi,gr,gi,pi, 
step, sumr, sumi pc, 
aintr,ainti 
J,n 

call «tep/ind(xi,step,J) 
sumr = O.ODO 
sumi = O.ODO 
x = O.ODO 

call rjrand(x,xi,qr,qi,aintr,aint pi) 
do 10 n=l,J 
x = x + step 

call rjrand(x,xi,qr,qi,bintr,binti,pi) 
sumr = sumr 4 aintr 4 4D0*bintr 
sumi = sumi 4 ainti 4 4DG*binti 
x — x + step 

cab r yran (/(x,xi,qr,qi, aintr, ainti, pi) 
sumr ~ sumr + airitr 
sumi ~~ sumi -f ainti 
10 continue 


double 

+ 

+ 

integer 


4 
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gr = step*sumr/3D0 
gi = step*sumi/3D0 

return 

end 


subroutine ryrani(x,xi,qr,qi,intr,inti,pi) 
double precision x,xi,qr,qi,intr,inti,pi, 


+ 

y,vr,vi,dvr,dvj,wlr,wli, 

+ 

dwlr,dwli,w2r,w2i,dw2r,dw2i 

+ 

c ,s , ar ,ai ,br ,bi , denomr, denomi , 

+ 

denomsq , recipr , recipi ,cr , ci . 

+ 

mult,int2r,int2i,intlr ) intli 


y = 0D0 
c = 0.5D0 

s = 0.5D0*dsgr£(3D0) 

call Fc?c*(x,y,vr,vi,dvr,dvi,wlr,wli,dwlr,dwli, 
+ w2r,w2i ) dw2r,dw2i) 

call cprod(-c,s,qr,qi,ar ; ai) 
call cprod(ar,ai,w2r,w2i,br,bi) 
denomr = dw2r - br 
denomi = dw2i - bi 
denomsq = denomr**2 + denomi**2 
recipr = denomr/ denomsq 
recipi = - denomi /denomsq 
ar = <fcos(x*xi/2D0) 
ai — - dstn(x*xi/2D0) 
call cprod (ar,ai, recipr, recipi, cr,ci) 
mult ~ dczp(-x*x i*s)/ ds^ri(pi) 
int2r — mult*cr 
int2i - mult*ci 

call cprod(qr ,qi t wlr ,v/li y br ,bi) 

denomr — dwlr - br 

denomi = dwli - bi 

denomsq = denomr**2 + denomi**2 

recipr = denomr/ denomsq 

recipi = - denomi /denomsq 

ar = <fcoa(x*xi) 

ai = dsin(x*x\) 

call cprtfd( ar,ai, recipr, recipi, cr,ci) 
mult = IDO / </s?rf(pi) 
intlr = mult*cr 
intli — rnult*ci 
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intr •= intlr + int2r 
inti = intli + int2i 


return 

end 


subroutine atepfind{x i,step,J) 

double precision xi,stepjay 
integer J 

if (<fo6s(xi) .It. 2D0) 

-I- then 

step = 0.05D0 

else 

step = 0.1D0/ daba(x i) 
end if 

jay = 3.77D0/ step 
J - jay 

return 

end 


subroutine ridgput(x i,qr.qi) 

double precision xi,qr,qi 

write (*,*) ’Program Ridgchek’ 
write (*,*) ’Version of February 1986’ 
write (*,*) 
write (*,*) 

write (*,*) ’What is magnitude of argument xi? ’ 
read (*,*) xi 

write (*,*) ’What is real part of q? ’ 
read (*,*) qr 

write (*,*) ’What is imaginary part of q? ’ 
read (*,*) qi 

return 

end 


subroutine prtndye(xi,qr,qi,gr,gi,g,lo88) 
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double precision xi,qr,qi,gr,gi,gJ[oss 


write (*,*) 
write (*,*) 
write (*,1) ’ xi = ’, xi, 

/i av a n < ^ 


vrite (*,1J ' xi = \ xi, * qr = qr, ’ qi = qi 
format (10X, A, F10.4, 10X, A, F10.4, A, “ 
vrite (*,*) 


write ^ , j 

write (*,3) ’ gr\ ’ gi’, ’ g’, ’ ioss’ 
format (12X,A,11X,A,10X,A,10X,A) 
write (*,*) 


F10.4) 


V 5 j 

write (*,101) gr,gi,g,loss 
write (*,*) 


^ , J 

101 format (lx, 4D16.8) 


return 

end 


Complex roots for creeping wave transcendental equation 


The program listed here evaluates the roots of the transcendental equation in Eq. (4). The 

present verson is temporary and allows the user to iterate refinements based on Newton’s method 
at the keyboard. 


program crecpwve 


Allan D. Pierce 
2/ 13/86 


double precision xrstart,xistart,qr,qi, 
+ xrfin,xifin,gr,gi 


call creeput (xrstart, xistart, qr,qi) 
call neu;tcrp(xrstartpcistart,qr ) qi,gr,gi) 
xrfin = xrstart - gr 
xifin = xistart - gi 

call prtcreep (xrstart, xistart, xrfin, xifin) 


stop 

end 


subroutine newtcrp(x, y,qr,qi,gi,gi) 


+ 

+ 

+• 


double preci ion x.y.qr.qi.gr.gi.ar.ai, 


vr,vi,dvr,dvi,wlr,wli, 
dwlr ) dwli,w 2 r,w 2 i,dw 2 r,dw 2 i, 
br,bi,denomr,denomi, 
denomsq ,recipr .recipi , 
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+ numr^mmi,3,c,er,ei,cr,ci 

s =0.5D0**9rt(3D0) 
c =0.5D0 

call cprod(x,y,-c ) s ) er,ei) 

call Fock (er,ei,vr,vi,dvr,dvi,wlr,wli,dw lr ,dw li, 
+ w2r,w2i,dw2r,dw2i) 

call cpro<i(qr,qi,c^,cr ) ci) 
call cprod( cr.ci.dvr.dvi.ar.ai) 
call cprod(er,ei,vr,vi,br,bi) 
denomr = br + ar 
denomi = bi + ai 

denomsq = denomr* *2 + denomi* *2 

recipr = denomr/ denomsq 

recipi ~ - denomi/ denomsq 

call c prod (cr, ci, vr,vi,ar,ai) 

numr = dvr -(• ar 

numi = dvi + ai 

call cprod(numr,numi,recipr,recipi ) er,ei) 
call cprod(er ) ei,-c,-3,gr,gi) 
return 
end 


4 



subroutine creeput(xr,xi,qr,qi) 

double precision xr,xi,qr,qi 

write (*,*) ’Program Creepwve’ 
write (*,*) ’Version of February 1986’ 
write (*,*) 

write (*,*) ’What is real part of q? 
read (*,*) qr 

write (*,*) ’What is imaginary part of q? 
read (*,*) qi 

write (*,*) ’What is real part of initial guess for root? 
read (*,*) xr 

write (*,*) ’What is imaginary part of initial guess? 
read (*,*) xi 

return 

end 

subroutine prtcreep(xrstart,xistart,xrfin,xifin) 
double precision xrstart,xistart,xrfin,xifin 
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